Pronóstico de OT (ETT) sin estacionalidad — 100 pasos sin regresores futuros

Objetivo: Seleccionar un dataset multivariante no financiero (ETT), comprobar que la variable objetivo OT (Oil Temperature) no tiene estacionalidad y construir ~5 modelos para predecir 100 pasos en el futuro sin acceso a regresores futuros desde el instante de predicción. Evaluar con KPIs (MAE, RMSE, MAPE), guardar el mejor modelo y mostrar una demo de uso.

Dataset: ETT-small ETTh1.csv (horario), columnas: date, HUFL, HULL, MUFL, MULL, LUFL, LULL, OT.

Fuente: https://github.com/zhouhaoyi/ETDataset

1. Preparación del entorno

import os
import sys
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from ydata_profiling import ProfileReport


from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor

from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.tsa.seasonal import STL
from statsmodels.tsa.api import ExponentialSmoothing
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf

import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers

from math import sqrt
from sklearn.metrics import mean_absolute_error, mean_squared_error

import joblib

# Local utils
sys.path.append(str(Path.cwd() / "src"))


BASE_DIR = Path('.').resolve().parent if (Path('.').resolve().name == 'notebooks') else Path('.').resolve()
DATA_RAW = BASE_DIR / 'data' / 'raw'
ARTIFACTS = BASE_DIR / 'artifacts'
REPORTS = BASE_DIR / 'reports'
FIGURES = REPORTS / 'figures'
RANDOM_STATE = 42

# Crear directorios si no existen
for p in [DATA_RAW, ARTIFACTS, REPORTS, FIGURES]:
    p.mkdir(parents=True, exist_ok=True)

2. Descarga o carga del dataset

import urllib.request

url = "https://raw.githubusercontent.com/zhouhaoyi/ETDataset/main/ETT-small/ETTh1.csv"
csv_path = DATA_RAW / "ETTh1.csv"
if not csv_path.exists():
    print("Descargando ETTh1.csv...")
    urllib.request.urlretrieve(url, csv_path)

df = pd.read_csv(csv_path, parse_dates=["date"])
df = df.sort_values("date").reset_index(drop=True)
df.head()

	date	HUFL	HULL	MUFL	MULL	LUFL	LULL	OT
0	2016-07-01 00:00:00	5.827	2.009	1.599	0.462	4.203	1.340	30.531000
1	2016-07-01 01:00:00	5.693	2.076	1.492	0.426	4.142	1.371	27.787001
2	2016-07-01 02:00:00	5.157	1.741	1.279	0.355	3.777	1.218	27.787001
3	2016-07-01 03:00:00	5.090	1.942	1.279	0.391	3.807	1.279	25.044001
4	2016-07-01 04:00:00	5.358	1.942	1.492	0.462	3.868	1.279	21.948000

3. Exploración y selección de la variable objetivo

df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 17420 entries, 0 to 17419
Data columns (total 8 columns):
 #   Column  Non-Null Count  Dtype         
---  ------  --------------  -----         
 0   date    17420 non-null  datetime64[ns]
 1   HUFL    17420 non-null  float64       
 2   HULL    17420 non-null  float64       
 3   MUFL    17420 non-null  float64       
 4   MULL    17420 non-null  float64       
 5   LUFL    17420 non-null  float64       
 6   LULL    17420 non-null  float64       
 7   OT      17420 non-null  float64       
dtypes: datetime64[ns](1), float64(7)
memory usage: 1.1 MB

No hay valores nulos, veamos un análisis más exhaustivo con ProfileReport

report = ProfileReport(df)
report

Summarize dataset:   0%|          | 0/5 [00:00<?, ?it/s]

100%|██████████| 8/8 [00:00<00:00, 48.87it/s]

Generate report structure:   0%|          | 0/1 [00:00<?, ?it/s]

Render HTML:   0%|          | 0/1 [00:00<?, ?it/s]

Podemos remarcar varias cosas sobre el análisis realizado con Profile Report.

En primer lugar, si observamos la matriz de correlación, las variables HUFL (High UseFul Load) y MUFL (Middle UseFul Load) están altamente correlacionadas. Lo mismo ocurre con HULL (High UseLess Load) y MULL (Middle UseLess Load). Por este motivo, eliminaremos las variables MUFL y MULL y nos quedaremos solamente con variables de alta carga y baja carga.

Por otro lado, también tenemos presencia de outliers que impiden una distribución normal.

df_filtered = df.copy()
df_filtered.drop(['MUFL']  ,axis=1, inplace=True, errors = 'ignore')
df_filtered.drop(['MULL'] ,axis=1, inplace=True, errors = 'ignore')

# Función para eliminar outliers
def remove_outliers(df: pd.DataFrame, columns: list, threshold: int =1.5) -> pd.DataFrame:
    """Elimina outliers de un DataFrame utilizando el método del rango intercuartílico (IQR).

    Args:
        df (pd.DataFrame): DataFrame del cual se eliminarán los outliers.
        columns (list): Lista de columnas a evaluar para la eliminación de outliers.
        threshold (int, optional): Umbral por el cual se modifica la sensibilidad del rango (a mayor, menos datos se eliminan). Defaults to 1.5.

    Returns:
        pd.DataFrame: DataFrame sin outliers.
    """
    for col in columns:
        Q1 = df[col].quantile(0.25)
        Q3 = df[col].quantile(0.75)
        IQR = Q3 - Q1
        lower_bound = Q1 - threshold * IQR
        upper_bound = Q3 + threshold * IQR
        df = df[(df[col] >= lower_bound) & (df[col] <= upper_bound)]
    return df

# Remove outliers from the filtered dataset
columns_to_check = ['HUFL', 'HULL', 'LUFL', 'LULL', 'OT']
df_filtered = remove_outliers(df_filtered, columns_to_check, 2)

Veamos de nuevo el report

report_filtered = ProfileReport(df_filtered)
report_filtered

Summarize dataset:   0%|          | 0/5 [00:00<?, ?it/s]

100%|██████████| 6/6 [00:00<00:00, 226.89it/s]

Generate report structure:   0%|          | 0/1 [00:00<?, ?it/s]

Render HTML:   0%|          | 0/1 [00:00<?, ?it/s]

Ahora ya no hay tanta presencia de outliers y siguen distribuciones gaussianas.

Representaciones interesantes para el estudio:

- Evolución de **OT** total.
- Evolución de **OT** por días en el mes.
- Evolución de **OT** por días en la semana.
- Evolución de **OT** por meses.

# Evolución de OT total
plt.figure(figsize=(12, 6))
plt.plot(df_filtered["date"], df_filtered["OT"], label="OT", color="blue", linewidth=1.5)
plt.title("Evolución de OT por Fecha", fontsize=16)
plt.xlabel("Fecha", fontsize=12)
plt.ylabel("Oil Temperature (OT)", fontsize=12)
plt.grid(alpha=0.3)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig(FIGURES / "ot_evolution.png", dpi=160)
plt.show()

Se observan períodos de mayor temperatura en los meses de verano, con patrones a corto tiempo. Sin embargo, se puede apreciar una tendencia negativa. Veamos más profundamente cada caso.

Plot de la media de OT por día por cada mes durante los dos años.

import math

# Extraemos mes-año, día de la semana y días completos sin hora
df_filtered['YearMonth'] = df_filtered['date'].dt.to_period('M')
df_filtered['Weekday'] = df_filtered['date'].dt.day_name()  # lunes, martes...
df_filtered['Day'] = df_filtered['date'].dt.date
# Si prefieres orden numérico de lunes=0,... domingo=6:
df_filtered['WeekdayNum'] = df_filtered['date'].dt.weekday

# Agrupamos por mes y día de la semana
monthly_daily_mean = df_filtered.groupby(['YearMonth', 'WeekdayNum'])['OT'].mean().reset_index()

# Agrupación por mes
monthly_mean = df_filtered.groupby('YearMonth')['OT'].mean().reset_index()
monthly_mean['YearMonth'] = monthly_mean['YearMonth'].dt.to_timestamp()

# Agrupación por día (pd.Series)
series_daily_mean = df_filtered.set_index('date').resample('D')['OT'].mean()
daily_mean = df_filtered.groupby('Day')[['OT', 'HUFL', 'HULL', 'LUFL', 'LULL']].mean().reset_index()

# Mapeamos el número a nombre del día para el plot ordenado
weekday_map = {0:'Mon', 1:'Tue', 2:'Wed', 3:'Thu', 4:'Fri', 5:'Sat', 6:'Sun'}
monthly_daily_mean['Weekday'] = monthly_daily_mean['WeekdayNum'].map(weekday_map)

# Lista de meses únicos ordenados
unique_months = sorted(monthly_daily_mean['YearMonth'].unique())

# Configuración de la grid de subplots
n_months = len(unique_months)
ncols = 4   # columnas en la grid (puedes ajustar)
nrows = math.ceil(n_months / ncols)

fig, axes = plt.subplots(nrows, ncols, figsize=(4*ncols, 3*nrows), sharey=True)

# Flatten para iterar fácilmente
axes = axes.flatten()

for i, ym in enumerate(unique_months):
    subset = monthly_daily_mean[monthly_daily_mean['YearMonth'] == ym]
    axes[i].bar(subset['Weekday'], subset['OT'])
    axes[i].set_title(str(ym))
    axes[i].set_xlabel("Weekday")
    if i % ncols == 0:  # solo la primera columna lleva ylabel
        axes[i].set_ylabel("Mean OT")

# Eliminar subplots vacíos
for j in range(i+1, len(axes)):
    fig.delaxes(axes[j])

plt.tight_layout()
plt.savefig(FIGURES / "ot_weekly_mean_per_month.png", dpi=160)
plt.show()

Gráficos de la temperatura del transformador durante cada día de cada mes

months = ['January','February','March', 'April', 'May','June', 
              'July', 'August', 'September', 'October', 'November', 'December']

for year in series_daily_mean.index.year.unique():
    subset = series_daily_mean[series_daily_mean.index.year == year]
    plt.figure(figsize=(12,4))
    for month in subset.index.month.unique():
        plt.plot(
            subset[subset.index.month == month].index.day,
            subset[subset.index.month == month].values,
            marker='o',
            label=str(months[month-1])
        )

    plt.title(f"Temperatura media diaria - {year}")
    plt.xlabel("Día del mes")
    plt.ylabel("OT (media diaria)")
    plt.legend()
    plt.grid(True)
    plt.savefig(FIGURES / f"mean_daily_ot_{year}", dpi=160)
    plt.show()

Al analizar la media de temperatura por día y por mes, se observa un incremento de la variable OT durante los meses de verano. Asimismo, al comparar periodos equivalentes en distintos años (por ejemplo, julio de 2016 frente a julio de 2017), se aprecia un comportamiento con patrones similares. No obstante, el aspecto más evidente de los gráficos es la tendencia descendente de la temperatura del aceite a lo largo del tiempo.

Este comportamiento constituye un primer indicio de que la variable objetivo OT no es estacionaria. Sin embargo, dado que únicamente se dispone de dos años de datos con frecuencia mensual, no resulta posible confirmar formalmente la no estacionariedad a partir de un patrón anual, pues sería necesario contar con una serie temporal más extensa. En este sentido, una alternativa más adecuada consiste en agrupar los datos por media diaria y comprobar empíricamente la no estacionariedad antes de avanzar en el modelado y la predicción de la temperatura diaria.

# Plot de media de temperatura por mes, de nuevo se puede observar la tendencia negativa.
plt.figure(figsize=(6,3))
plt.plot(monthly_mean['YearMonth'], monthly_mean['OT'], marker='o')
plt.title("Media de temperatura por mes (Jul 2016 - Jul 2018)")
plt.xlabel("Mes")
plt.ylabel("Mean OT")
plt.grid(True)
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

4. Comprobación de no estacionariedad de OT

Comprobación estadística de que OT es no estacionario

y = daily_mean["OT"].values


# STL (periodo de 30 dias si fuera estacional mensual)
stl = STL(daily_mean["OT"], period=30)  # probar periodo mensual teórico
res = stl.fit()
fig = res.plot()
fig.suptitle("Descomposición STL (periodo tentativo=30)")
plt.tight_layout()
plt.savefig(FIGURES / "STL.png", dpi=160)
plt.show()

# ACF para buscar picos periódicos (30, 60, 180, etc.)
fig = plt.figure()
plot_acf(y, lags=180)
plt.title("ACF de OT (hasta lag 24)")
plt.tight_layout()
plt.savefig(FIGURES / "ACF.png", dpi=160)
plt.show()

# Pruebas ADF y KPSS
adf_stat, adf_p, *_ = adfuller(y, autolag="AIC", regression="ct")
kpss_stat, kpss_p, *_ = kpss(y, regression="ct", nlags="auto")

print(f"ADF p-value: {adf_p:.4f}  (H0: raíz unitaria)")
print(f"KPSS p-value: {kpss_p:.4f} (H0: estacionariedad)")

# Comentario automático
if adf_p < 0.05:
    print("ADF sugiere rechazo de H0: no hay raíz unitaria (posible estacionaria).")
else:
    print("ADF no rechaza H0: posible no estacionaria.")
if kpss_p < 0.05:
    print("KPSS rechaza H0: posible no estacionaria (tendencia).")
else:
    print("KPSS no rechaza H0: estacionaria.")

<Figure size 640x480 with 0 Axes>

ADF p-value: 0.7302  (H0: raíz unitaria)
KPSS p-value: 0.0100 (H0: estacionariedad)
ADF no rechaza H0: posible no estacionaria.
KPSS rechaza H0: posible no estacionaria (tendencia).

C:\Users\belop\AppData\Local\Temp\ipykernel_2876\2492847394.py:23: InterpolationWarning: The test statistic is outside of the range of p-values available in the
look-up table. The actual p-value is smaller than the p-value returned.

  kpss_stat, kpss_p, *_ = kpss(y, regression="ct", nlags="auto")

Ambos test corroboran la no estacionariedad de la variable OT

5. Partición temporal y normalización

Uso de daily_mean, datos agrupados por media de horas en cada día

# Reservamos los últimos 100 puntos para evaluación de predicción multi-step
import joblib


TEST_HORIZON = 100

# División de train y test
train = daily_mean.iloc[:-TEST_HORIZON].copy()
test = daily_mean.iloc[-TEST_HORIZON:].copy()

y_train = train["OT"].values.astype(float)
y_test = test["OT"].values.astype(float)

# Escalado (solo con train)
scaler = StandardScaler()
y_train_sc = scaler.fit_transform(y_train.reshape(-1,1)).ravel()
y_test_sc = scaler.transform(y_test.reshape(-1,1)).ravel()

# Guardado de scaler
joblib.dump(scaler, ARTIFACTS / 'scaler_ot_daily_mean.pkl')

print(len(y_train), len(y_test))

602 100

6. Estrategia de predicción 100 pasos (100 días)

# Usaremos una estrategia recursiva: entrenamos modelos one-step-ahead con lags de OT.
# Después, a partir del último histórico conocido, predecimos el siguiente valor y lo alimentaremos
# (feed-forward) como si fuera observado para predecir el siguiente, y así sucesivamente.

N_LAGS = 302  # 2 días como ventana inicial (ajustable)
Xtr, ytr = [], []
for i in range(N_LAGS, len(y_train_sc)):
    Xtr.append(y_train_sc[i-N_LAGS:i])
    ytr.append(y_train_sc[i])
Xtr = np.array(Xtr)
ytr = np.array(ytr)

print(Xtr.shape, ytr.shape)

(300, 302) (300,)

7. Modelos a comparar

from typing import Dict, Callable
models: Dict[str, object] = {}

# 7.1 Naive (baseline): persiste último valor
class NaiveLast:
    def __init__(self, n_lags=48):
        self.n_lags = n_lags
    def fit(self, X, y):
        return self
    def predict(self, X):
        # devuelve el último valor de la ventana (identidad)
        return X[:, -1]
models["NaiveLast"] = NaiveLast(n_lags=N_LAGS)

# 7.2 ARIMA (entrenado sobre y_train_sc). Elegimos pequeño grid.
# Para 1-step, predecimos recursivamente fuera de muestra.
arima_order_candidates = [(2,0,2), (3,0,2), (2,1,2), (3,1,2)]
best_aic = np.inf
best_order = None
best_arima = None
for order in arima_order_candidates:
    try:
        m = ARIMA(y_train_sc, order=order).fit(method_kwargs={"warn_convergence":False})
        if m.aic < best_aic:
            best_aic = m.aic
            best_order = order
            best_arima = m
    except Exception as e:
        pass
class ARIMAWrapper:
    def __init__(self, model_fit, n_lags=48):
        self.model_fit = model_fit
        self.n_lags = n_lags
    def fit(self, X, y):
        return self
    def predict(self, X):
        # 1-step ahead pred usando último estado del modelo entrenado: proxy simple
        # En fase recursiva fuera del set, usaremos simulate via rolling
        # Aquí devolvemos el último residuo + forecast 1-step de entrenamiento como aproximación
        return np.full((X.shape[0],), self.model_fit.forecast(1)[0])
models[f"ARIMA{best_order}"] = ARIMAWrapper(best_arima, n_lags=N_LAGS)

# 7.3 Holt (Exponential Smoothing)
class HoltWrapper:
    def __init__(self, y_train_sc, n_lags=48):
        self.model = ExponentialSmoothing(y_train_sc, trend='add', seasonal=None).fit()
        self.n_lags = n_lags
    def fit(self, X, y):
        return self
    def predict(self, X):
        # usar forecast 1-step constante como proxy para lotes
        f = self.model.forecast(1)[0]
        return np.full((X.shape[0],), f)
models["Holt"] = HoltWrapper(y_train_sc, n_lags=N_LAGS)

# 7.4 RandomForest
rf = RandomForestRegressor(n_estimators=300, max_depth=12, random_state=42, n_jobs=-1)
models["RandomForest"] = rf

# 7.5 XGBoost
xgb = XGBRegressor(
    n_estimators=600,
    max_depth=8,
    learning_rate=0.05,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42,
    n_jobs=-1,
    reg_lambda=1.0,
)
models["XGBoost"] = xgb

# 7.6 LSTM (opcional; puede ser más lento)
try:

    def build_lstm(n_lags=48):
        model = keras.Sequential([
            layers.Input((n_lags, 1)),
            layers.LSTM(64, return_sequences=False),
            layers.Dense(32, activation='relu'),
            layers.Dense(1)
        ])
        model.compile(optimizer=keras.optimizers.Adam(1e-3), loss='mse')
        return model
    models["LSTM"] = "placeholder"
except Exception as e:
    print("TensorFlow no disponible:", e)

list(models.keys())

c:\Users\belop\Documents\Trabajo\Capgemini\CBL-TestDataScience-2\.prueba_env\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:966: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
  warn('Non-stationary starting autoregressive parameters'
c:\Users\belop\Documents\Trabajo\Capgemini\CBL-TestDataScience-2\.prueba_env\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.
  warn('Non-invertible starting MA parameters found.'

['NaiveLast', 'ARIMA(2, 1, 2)', 'Holt', 'RandomForest', 'XGBoost', 'LSTM']

8. Entrenamiento (one-step)

results = {}

# Preparación para LSTM
Xtr_lstm = Xtr.reshape((Xtr.shape[0], Xtr.shape[1], 1))

for name, model in models.items():
    if name == "LSTM":
        m = build_lstm(N_LAGS)
        m.fit(Xtr_lstm, ytr, epochs=10, batch_size=128, verbose=0, validation_split=0.1)
        models[name] = m  # replace placeholder
        results[name] = {"status": "trained"}
    else:
        model.fit(Xtr, ytr)
        # store attribute n_lags for demo script
        if hasattr(model, "n_estimators") or hasattr(model, "predict"):
            setattr(model, "n_lags", N_LAGS)
        results[name] = {"status": "trained"}

results

{'NaiveLast': {'status': 'trained'},
 'ARIMA(2, 1, 2)': {'status': 'trained'},
 'Holt': {'status': 'trained'},
 'RandomForest': {'status': 'trained'},
 'XGBoost': {'status': 'trained'},
 'LSTM': {'status': 'trained'}}

9. Predicción recursiva de 100 pasos

def recursive_forecast(model, last_hist_scaled, n_steps=100, is_lstm=False):
    hist = list(last_hist_scaled)  # last N_LAGS values, scaled
    preds = []
    for _ in range(n_steps):
        x = np.array(hist[-N_LAGS:])
        if is_lstm:
            x_in = x.reshape(1, N_LAGS, 1)
            yhat = float(model.predict(x_in, verbose=0).ravel()[0])
        else:
            x_in = x.reshape(1, -1)
            yhat = float(model.predict(x_in).ravel()[0])
        preds.append(yhat)
        hist.append(yhat)
    return np.array(preds)

hist_window = y_train_sc[-N_LAGS:]

preds_dict = {}
for name, model in models.items():
    if name.startswith("ARIMA"):
        # For ARIMA/Holt we use their 1-step forecast in a recursive fashion by refitting rolling.
        # Simple strategy: persist last value or use a small wrapper; we fallback to Naive-like behavior
        # if direct recursive with internal state is not straightforward in batch.
        # Here we approximate with recursive using last window + model 1-step.
        preds_sc = recursive_forecast(model, hist_window, n_steps=100, is_lstm=False)
    elif name == "Holt":
        preds_sc = recursive_forecast(model, hist_window, n_steps=100, is_lstm=False)
    elif name == "LSTM":
        preds_sc = recursive_forecast(model, hist_window, n_steps=100, is_lstm=True)
    else:
        preds_sc = recursive_forecast(model, hist_window, n_steps=100, is_lstm=False)
    preds = scaler.inverse_transform(preds_sc.reshape(-1,1)).ravel()
    preds_dict[name] = preds

# Añadimos baseline Naive directo usando último valor observado
last_obs = y_train[-1]
preds_dict["NaiveRepeat"] = np.repeat(last_obs, 100)

list(preds_dict.keys())

['NaiveLast',
 'ARIMA(2, 1, 2)',
 'Holt',
 'RandomForest',
 'XGBoost',
 'LSTM',
 'NaiveRepeat']

Predicción directa (no recursiva)

# ARIMA: re-fit sobre todo el y_train_sc si no lo has guardado
arima_fit = ARIMA(y_train_sc, order=best_order).fit()
arima_preds_sc = arima_fit.forecast(steps=TEST_HORIZON)  # devuelve array (h,)
arima_preds = scaler.inverse_transform(arima_preds_sc.reshape(-1,1)).ravel()

# Holt (ExponentialSmoothing)
holt_model = ExponentialSmoothing(y_train_sc, trend='add', seasonal=None).fit()
holt_preds_sc = holt_model.forecast(TEST_HORIZON)
holt_preds = scaler.inverse_transform(holt_preds_sc.reshape(-1,1)).ravel()

preds_dict["ARIMA_no_Recursive"] = arima_preds
preds_dict["Holt_no_Recursive"] = holt_preds

10. Evaluación con KPIs y comparación

def compute_metrics(y_true, y_pred):
    mae = mean_absolute_error(y_true, y_pred)
    rmse = sqrt(mean_squared_error(y_true, y_pred))
    mape = np.mean(np.abs((y_true - y_pred) / np.maximum(np.abs(y_true), 1e-8))) * 100.0
    return mae, rmse, mape

scores = []
for name, preds in preds_dict.items():
    mae, rmse, mape = compute_metrics(y_test, preds)
    scores.append((name, mae, rmse, mape))

scores_df = pd.DataFrame(scores, columns=["Modelo", "MAE", "RMSE", "MAPE"]).sort_values("RMSE")
scores_df.to_csv(REPORTS / 'scores_df.csv', index=False)
scores_df

	Modelo	MAE	RMSE	MAPE
4	XGBoost	1.992506	2.504265	25.776200
5	LSTM	2.481418	2.942341	30.162100
1	ARIMA(2, 1, 2)	3.091598	3.592667	34.826173
7	ARIMA_no_Recursive	3.240885	3.754269	35.997380
0	NaiveLast	3.376388	3.901175	37.013833
6	NaiveRepeat	3.376388	3.901175	37.013833
2	Holt	3.395377	3.920474	37.174691
3	RandomForest	3.514913	4.056638	37.826475
8	Holt_no_Recursive	4.545434	5.090525	48.913669

Gráfico comparativo (real vs. mejor modelo)

best_name = scores_df.iloc[0]["Modelo"]
best_preds = preds_dict[best_name]

plt.figure(figsize=(12,6))
plt.plot(test["Day"], y_test, label="Real")
plt.plot(test["Day"], best_preds, label=f"Predicción - {best_name}")
plt.title(f"Pronóstico 100 pasos — Mejor modelo: {best_name}")
plt.xlabel("Fecha")
plt.ylabel("OT")
plt.legend()
plt.tight_layout()
plt.savefig(FIGURES / "best_forecast.png", dpi=160)
plt.show()

best_name

'XGBoost'

11. Conclusión

Las series temporales no estacionarias suelen presentar una mayor complejidad en el proceso de predicción de evolución futura. A diferencia de las series estacionarias, en las que existen patrones claros y recurrentes, las series no estacionarias muestran una elevada irregularidad en los valores, lo que puede conllevar cambios abruptos e impredecibles.

En el caso del presente análisis, el dataset correspondiente a dos transformadores eléctricos —con atributos asociados a carga útil/inútil y carga alta/baja— contiene una variable objetivo (Oil Temperature) que se comporta de manera claramente no estacionaria tanto en periodos mensuales como, aparentemente, en escalas anuales.

Para abordar la tarea de predicción de los 100 periodos temporales establecidos, se llevó a cabo un proceso de visualización y limpieza de los datos. Posteriormente, a partir de un análisis exploratorio mediante gráficos y estadísticas descriptivas, se optó por realizar la predicción de la temperatura media diaria, considerando los últimos 100 días como conjunto de prueba.

Se evaluaron siete algoritmos de predicción: NaiveLast, ARIMA(2,1,2), Holt, Random Forest, XGBoost, LSTM y NaiveRepeat. En el caso de Holt y ARIMA, se implementaron dos configuraciones: predicciones recursivas de un día en un día y predicción directa de los 100 días. El tamaño de ventana empleado para el entrenamiento fue de 300 días, con el fin de aprovechar al máximo la limitada cantidad de datos disponible.

Entre todos los modelos, XGBoost obtuvo los mejores resultados en términos de métricas de evaluación (MAE, MAPE y RMSE). Sin embargo, a pesar de superar al resto de algoritmos, su rendimiento no alcanza niveles óptimos: el MAPE obtenido fue del 25%, cuando lo recomendable en tareas de predicción suele ser un valor inferior al 5%.

No obstante, el análisis visual de las predicciones revela que XGBoost consigue capturar la tendencia global de la serie, reproduciendo adecuadamente los incrementos de temperatura en periodos específicos. Un ejemplo representativo se observa alrededor del 10 de abril de 2018, donde el modelo refleja el aumento de la temperatura del aceite atribuido a la proximidad del verano, lo que indica una capacidad para identificar patrones estacionales de carácter aproximado.

En conclusión, aunque el modelo de predicción basado en XGBoost muestra utilidad para anticipar la tendencia media mensual de la temperatura, las fuertes fluctuaciones de la serie limitan su capacidad de estimar con precisión los valores diarios.

11. Guardado del mejor modelo y demostración de carga

import joblib

best_model_obj = models.get(best_name, None)
# Guardamos el mejor modelo + scaler juntos
bundle_path = ARTIFACTS / "best_model.joblib"
joblib.dump({"model": best_model_obj, "scaler": scaler}, bundle_path)
bundle_path

WindowsPath('C:/Users/belop/Documents/Trabajo/Capgemini/CBL-TestDataScience-2/artifacts/best_model.joblib')

Puedes ejecutar la demo de inferencia desde terminal:

python scripts/demo_inferencia.py